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Task  (i):  Collaborative  development  of  a  circuit  system  that  is  capable  of  accepting  an  input  pulse  and 
converting  it  to  the  traveling  solitary  wave  pulse  seen  in  a  granular  alignment  ( Solitary  Wave  Chip 
Problem). 


The  write  up  below  describes  our  accomplishments  in  constructing  a  Simulink  version  of  the  granular 
chain. 


Simulink  Modeling  of  Granular  Chains  (Robert  W.  Newcomb,  Department  of  Electical  and  Computer 
Engineering,  University  of  Maryland,  College  Park,  Maryland  20742  USA  (newcomb@eng.umd.edu) 
Surajit  Sen,  Physics  Department,  State  Univesity  of  New  York,  Buffalo,  New  York,  14260 
(sen@buffalo.edu)) 

Abstract — After  a  review  of  the  coupled  Newton’s  equations  for  a  granular  alignment,  the  equations  are 
put  into  block  diagrams  of  Simulink.  Simulink  simulations  are  given  for  22  grain  systems  with  cubic  and 
for  Hertz  potential  energy.  The  expected  granular  solitary  waves  are  seen  in  the  simulations. 


1.  Introduction 


The  types  of  grains  under  consideration  comprise  a  one-dimensional  chain  of  symmetric  elastic  grains 
such  that  an  input  pulse  travels  through  compression  along  the  chain.  By  experiment  [1],  [2],  through 
simulations  [3],  and  by  series  approximations  [4],  [5],  the  pulses  are  known  to  be  able  to  form  into  solitary 
waves  and  since  action  potentials  are  solitary  waves,  these  are  similar  to  the  signals  used  by  biological 
neurons  [6,  p.  42]  and  of  considerable  interest  for  mimicking  neural  information  processing.  Therefore, 
these  grains  can  be  seen  as  an  alternate  means  of  forming  the  pulses  used  in  silicon  based  pulse  coded 
neural  networks  [7].  Alternatively,  their  equations  can  be  put  into  a  form  which  allows  for  an  equivalent 
transistor  structure  having  the  key  properties  of  the  elastic  spherical  grains.  Consequently,  with  an  ultimate 
goal  of  mimicking  the  grains  behavior  in  transistor  circuits,  in  this  paper  we  present  a  Simulink  model  of 
these  grains  in  a  form  that  allows  for  future  conversion  into  VLSI  circuits. 


Figurel .  Chain  of  grains  of  radii  R  placed  between  rigid  walls. 
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II.  Describing  Equations 

Figure  1  gives  a  one-dimensional  representation  of  the  granular  spheres  which  we  here  assume  all  have 
the  same  radius  R.  We  consider  N  grains  with  q  being  the  coordinate  of  the  center  of  the  ith  grain.  For  i=l 
an  external  impulse-like  force  is  assumed  applied  while  for  the  final  grain,  at  xN,  a  rigid  wall  is  assumed. 
We  use  the  Hamiltonian,  H(p,q)  representation  where  p  =  momentum  N-vector  and  q  =  position  N-vector 
and  H  is  the  sum  of  the  kinetic  and  the  potential,  V(.,.),  energies.  Thus 

H(p,q)  =  z  (-^p?  +  v(qi_i»qi»  (la)> 

0q. 

Pi=mlF’  (lb)> 

V(qi_1,qi)  =  k[(qi  l  +R)-(qi  -R)]^+1  (lc) 

Here  q  is  the  position  of  the  center  of  the  ith  grain  measured  from  an  origin  qrR  =  0.  The  potential 
energy  depends  on  the  overlap,  2R-(q1  -  qM),  of  two  adjacent  grains  if  positive  (and  is  zero  if  there  is  no 
overlap).  If  the  rest  position  is  qio  and  the  displaced  position  is  x,  =  q,  -  qio  then  the  overlap  is  x,.,-x„  which 
gives  the  potential  energy  V  if  positive  with  V  being  zero  if  there  is  no  overlap.  So,  following  [3]  the 
modified  symbol  [x]+  =  (x+|x|)/2=(l+sign(x))x/2  is  used  in  (1)  to  designate  x  if  x  >  0  and  zero  if  x  <  0.  The 
power  r+1  is  due  to  Hertz  [8]  and  known  to  be  5/2  for  elastic  spheres  though  we  simulate  with  other  r  as 
well  [9];  especially  r  =  2  is  convenient  for  analytic  investigations.  We  set  up  the  equations  with  r  as  a 
parameter  which  then  becomes  easy  to  change  in  Simulink.  The  mass  of  a  grain  is  m  and  k  comprises 
various  constants  including  Young’s  modulus. 

Simulionk  realizations  are  most  easily  obtained  through  the  state  variable  equations  which  in  this  case 
are  the  Hamilton  differential  equations. 


^USH=^d  (2a) 

<9t  dp.  mPi  m  dt  ^  ^ 


8R-  m^L 

5qi“mat2 


By  normalizations,  the  choice  of  x;,  and  introducing  possible  loss  (by  the  parameter  kloss)  we  recast 
these  into  the  following  state-variable  form  which  are  the  actual  ones  we  put  into  Simulink  in  the  following 
paragraphs. 
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(3a), 


dx. 

—L  =  z 
dt  i 

dx. 

dz.  d(— 77-) 

_ i_  _  dt 

dt  dt 


(3b) 


=  (-M  {((1  +  sign(x .  _  1  -  x.  ))(x.  _  1  -  x.  ))r 
2m 

-((l  +  sign(x. -x.  ,))(x. -x.  ,))r}-k,  z.. 

vv  6  v  1  i  +  l//v  1  i  +  l"  J  loss  1 

Equations  (3)  are  for  i  =  2,  N  while  at  i  =  1  an  additive  input  term,  f(t),  is  to  be  added  and  the  x0  term 
omitted  while  at  i  =  N  a  fixed  boundary  is  imposed  by  fixing  xN+1  =  (2N+1)R.  However,  since  the 
differences  of  position  hold,  the  x,  can  be  interpreted  as  the  incremental  displacement  of  the  center  of  the  ith 
sphere.  For  solitary  waves  of  velocity  c  we  have  x,(t)  =  u(x,-ct).  Following  normalizations  of  Chatterjee  [3], 
this  gives  the  second  order  differential  equation  for  the  solitary  wave 


u  =  [u(t+l)-u(t)£ -[u(t)-u(t-l)fj_.  (4) 

From  these  Chatterjee  shows  MatLab  simulations  indicating  the  existence  of  solitary  waves  while  Sen  & 
Manciu  [4]  give  a  series  solution  approximation.  In  detail,  [5],  with  a  =  xi-ct  and  n  a  parameter. 

u(a,  n))  =  ^  (1  -  tanh(F(a(n)))  (5a), 

F(«,n))  =  \  I  C2o+l(n)a2q+1  (5b)- 

2q=o  ^q+1 

The  C2q+i  have  been  evaluated  for  q  =  0,...,5  and  the  results  shown  to  be  solitary  type  waves  (see  Fig.  3.1 
of  [5]).  In  short  these  grains  are  known  to  support  solitary  waves. 

Consequently  we  know  that  we  can  obtain  solitary  waves  from  the  state  variable  equations  (3)  so  it  is  to 
them  we  turn  for  possible  transistor  realization.  Toward  that  we  obtain  next  a  suitable  block  diagram 
realization. 


Ill.  Simulink  Block  Diagrams 

Although  we  simulate  for  much  larger  N,  for  convenience  of  illustration  Fig.  2  shows  a  Simulink  block 
diagram  for  N=5  stages  of  grains  with  the  sub-blocks  for  i=2,3,4  being  given  in  detail  by  Fig.  3  where  we 
have  allowed  for  different  grain  radii  but  choose  R=0  for  equal  size  grains  and  x  as  displacement  around 
equilibrium  (the  choice  of  x;=qi  is  possible  with  nonzero  R  in  this  setup).  The  input,  i=l,  and  output, 
i=N=5,  stages  are  given  in  Figs.  4  &  5.  Figure  3  realizes  equations  (3)  while  the  input  and  output  stages, 
are  simple  modifications  reflecting  their  different  loading.  In  Fig.  2  an  input  pulse  is  applied  on  the  left  to 
the  input  stage. 
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Figure  2.  Simulink  5  grains  Simulink  block  connection 


Figure  3.  Simulink  ith  internal  grain  stage 
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Iipitiectbi 


Figure  5.  Detailed  Simulink  output  grain  stage,  i=N 

IV.  Simulation  Results 

Figure  6(at  the  end  of  the  paper)  gives  a  plot  of  dx/dt  at  the  fourth  and  the  21st  stages  for  r=2  and  m/k 
normalized  to  1,  showing  solitary  waves  as  well  as  their  reflection  from  the  N=22  end  wall.  For  Fig.  6  a 
square  input  pulse  of  amplitude  10'6  and  pulse  width  t=5xl0.  This  results  in  a  traveling  wave  of  amplitude 
3.5x10  9  with  a  delay  of  t=2,500  at  the  4th  grain,  for  r=2.  A  reflected  wave  can  be  seen  at  t=32,000  with  the 
solitary  pulse  arriving  at  the  21st  grain  at  t=  15,000,  all  in  normalized  time. 

V.  Discussion 

For  obtaining  Simulink  models  we  have  put  the  grains  differential  equations  into  state  variable 
form,  (3)  above.  From  these  we  are  able  to  set  up  block  diagrams  which  use  only  integrators,  multipliers, 
square  roots,  and  summers.  These  are  conveniently  put  into  Simulink  through  which  we  have  again 
shown  that  solitary  signals  can  be  generated.  In  the  case  of  grains  satisfying  the  Flertz  potentials  these 
blocks  necessitate  square  roots  in  obtaining  the  3/2  power.  Flowever,  from  the  simulations  we  obtain 
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similar  results  for  powers  of  r=2  as  well  as  3/2,  though  with  a  different  time  scaling  as  for  example  the  4th 
stage  peak  occurs  at  t=200  for  r  =  3/2.  We  have  normalized  m/k=l  and  considered  x,  as  the  change  in 
displacement,  but  we  can  consider  it  alternatively  as  the  absolute  center  position,  q,,  of  the  ith  grain;  we 
have  chosen  the  former  for  the  given  block  diagrams  (which  also  use  R=1  for  the  right  wall).  In  the 
Simulink  system  we  have  the  capability  of  Figure  6.  Simulation  result:for  r=2;  solid  line  dx/dt  at  4thgrain  ; 
dashed  line  dx/dt  at  21st  grain  both  for  N=22  grain  using  any  real  r  as  well  as  the  added  possibility  of 
including  loss,  though  it  is  not  present  in  the  basic  grains  equations.  The  above  can  be  generalized  to 
allow  for  different  radii  and  for  different  materials  of  different  individual  grains  within  the  system,  and 
this  is  allowed  for  by  separate  ports  in  the  individual  grain  cell  blocks. 

Some  useful  additional  references  are  included  [10-20],  For  example,  there  are  other  effects  which  can  be 
included,  one  of  which  uses  the  “coefficient  of  restitution”  [11],  while  generalization  to  higher 
dimensions  is  possible,  such  as  for  sand  at  the  beach.  Also  the  above  equations  are  normalized  though 
denormalization  is  easily  carried  out.  ([20]]  gives  material  constants  for  various  materials). 
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Task  (ii):  Developing  a  heterogeneous  granular  alignment  system  that  can  compress  or  dilate  propagating 
solitary  waves  ( Solitary  Wave  Width  Control  Problem). 


Width  of  the  Solitary  wave  in  granular  alignments  (Diankang  Sun,  Yoichi  Takato,  Nicholas  DeMeglio 
and  Surajit  Sen,  SUNY  Buffalo) 

Solitary  waves  (we  denote  solitary  wave  as  SW  throughout)  that  are  about  5  grain  diameters  wide 
naturally  form  in  a  chain  of  unloaded  elastic  spheres  that  repel  with  overlap  □  as 


V(SiM)  =  aSii+l" ,  8>  0,  where  for  spheres  n  =  5/2. [Nesterenko  1983,  Lazaridi  and  Nesterenko  1985, 


Sen  and  Manciu  2001,  Sokolow,  et  al,  2007,  Sen,  et  al,  2008]  For  n  — »  2  and  n  — »  oo,  a  one-sided 
harmonic  ID  chain  of  cylindrical  or  disk  shaped  beads  and  a  chain  of  particles  with  hard-sphere-like 
repulsion  are  approached  respectively.  It  is  known  that  for  n  >  2  an  impulse  propagates  in  these  systems 
as  solitary  waves  [Nesterenko  1983,  LN  1985,  Sen  et  al,  2008]  The  width  W  of  these  solitary  waves  are 
known  to  depend  only  upon  n  but  a  detailed  understanding  of  W(n  )  has  not  been  available.  Here  we  use 
numerical  and  analytical  methods  to  study  how  the  SW  width  depends  on  n.  In  the  numerical  study,  a 
geometric  tool  is  employed  to  calculate  the  total  time  averaged  kinetic  energy  of  the  SW  by  using  the 
virial  Theorem.  The  idea  is  to  construct  an  isosceles  triangle  area  such  that  its  area  equals  that  of  the  SW. 
Hence  the  base  of  the  triangle  gives  the  width  of  the  SW.  We  then  analytically  explore  the  width  of  the 
SW.  Here  we  present  the  results  of  a  detailed  computational  and  theoretical  study  to  propose, 

W{n)  cc  (n  -  2)  “  ,  where  a  ~  1/3  . 

Introduction 

Now  consider  an  alignment  of  elastic  grains  that  are  held  between  fixed,  rigid  walls  that  would  perfectly 
reflect  energy.  The  grain-grain  potential  is  described  for  convenience  as 


V(Su+l)  =  a„S.Mn,  S>  0, 


(1) 


9 


the  geometric  properties  of  the  grain,  Y  and  □  are  the  Young’s  modulus  and  Poisson’s  ratio,  and 
□  □  □  denotes  the  overlap  of  two  adjacent  grains. 

We  start  by  considering  the  equation  of  motion  of  a  bead  in  the  chain, 

m  =  nan  {  [A,  +  ui-\  (0  -  ut  (Of*  -  [A,  +  w,  (0  -  uM  ^  }  (2) 

where  d|  is  the  pre -compression  effected  on  the  grain  in  the  chain.  Sound  propagation  is  not  admissible 
when  there  is  no  preloading,  i.e.,  when  A;  =  0  .  This  is  because  when  n  >  2,  the  right  hand  side  (RHS)  of 

the  equation  cannot  be  Taylor  expanded.  So  there  is  no  quadratic  term  that  is  admissible  on  the  RHS, 
which  implies  that  there  are  no  harmonic  oscillations  [Nesterenko,  1983,  Nesterenko,  1985,  Chatterjee, 
Sen,  1999,  2001],  Hence  no  sound  propagation  is  possible.  ut  denotes  the  displacement  of  the  grain  i  from 
its  equilibrium  position.  High-precision  numerical  studies  show  that  the  discrete  chain  systems  admit 
robust  solitary  waves  propagating  through  the  chain  with  a  unique  wave  width.  These  studies  also  reveal 
that  the  wave  width  depends  only  on  n  (Manciu  et  al,  1999a).  A  perturbation  that  is  sufficiently  weak  such 
that  only  elastic  compressions  of  grains  are  involved  results  in  the  eventual  formation  of  a  solitary  wave 
(SW)  as  discussed  in  Sokolow  et  al  2007. 

Earlier  work  of  Nesterenko  has  claimed  using  continuum  theory  that 

w  =  .  (3) 

n-2  V  6 

So,  when  » 2,  W—> oo  (no  SW).  This  is  an  expected  result  because  n  =  2  is  the  acoustic  case.  When 

/7 — >oo,  W  ~  2 R,  which  is  the  extreme  hard  sphere  limit  when  W  shrinks  to  its  minimum. 

Earlier  numerical  study  using  exhaustive  dynamical  simulations  (see  Fig.  4  of  Manciu,  Sen  and  Hurd, 
1999)  suggest  that  Nesterenko’s  formula  for  W may  need  some  improvement. 
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The  correction  is  important  because  the  width  formula  is  sometimes  used  to  find  the  effective  “n”  in  the 
experiments.  Hence  this  formula  is  potentially  useful  in  designing  granular  meta-materials  (Daraio, 
private  communication). 

Numerical  study 

We  carried  out  our  simulations  by  setting  the  masses  of  the  grains  m  =  2.3 14  x  10”' kg,  n  =  set  of  values 
from  2.001  to  7.0,  potential  prefactor  a  =  52.442  N/mm1 5  for  all  «’ s.  Although  a„  depends  on  n  as  we 
have  seen  in  Sun,  Daraio  and  Sen,  Phys.  Rev.  E  83,  066605  (2011),  the  reason  we  can  use  one  constant  is 
because  each  run  is  independent  of  the  other  such  that  the  underlying  physics  does  not  change.  Further, 
since  we  expect  W  to  depend  on  “n”,  the  choice  of  the  value  of  a„  is  a  matter  of  convenience.  We  also  set 
B  =  0  so  there  are  no  harmonic  terms.  We  used  the  velocity  Verlet  algorithm  for  integrating  the  coupled 
equations  of  motion.  The  integration  time  step  was  taken  to  be  Dt  =  10"5  Ds.  In  the  simulation  a  delta 

1  2 

function  perturbation  with  Erot=—mv0  =0.115715/  was  initiated  at  grain  1  at  t=0.  We  set 

dissipation  to  zero.  Energy  was  constant  through  the  106  steps  of  our  simulations  to  one  part  in  106. 

We  observed  the  formation  of  a  single  solitary  wave  —  a  results  first  reported  by  Sokolow  et  al. 
(Sokolow  2007).  Snapshots  of  a  propagating  SW  at  two  different  times  are  shown  below  in  Fig.  1,  where 
KE  of  the  grains  vs.  grain  positions  are  shown  for  a  propagating  SW  at  two  different  times.  The  center 
grain  is  well  defined  when  the  wave  consists  of  an  odd  number  of  grains  whereas  the  center  is  roughly  at 
the  edges  of  the  two  center  grains  for  a  SW  sitting  on  an  even  number  of  grains.  Nevertheless,  the  total 
energy  does  not  disperse  even  though  the  SW  oscillates  between  the  two  widths.  To  choose  either  an  odd 
or  an  even  numbered  wave  width  to  be  the  solitary  wave  width  is  then  a  choice  of  definition.  See  Fig.  1. 
we  have  chosen  the  odd  numbered  solitary  wave  width  to  represent  the  width  here. 
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occillation 


OcciHation  depends  on  the 


ivai/e  traveling  direction 


nature  of  the  discrete  system. 


It  is  a  characteristic  feature 
of  the  solitary  wave  of  this 


V 


kind. 


$ 


odd  £  of  solitary  wave 


width 


evev  tt  of  solitary  irave  width 
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Fig.  1.  Solitary  wave  is  traveling  with  a  characteristic  fluctuation  in  width. 


Fig.  1  shows  that  the  width  of  the  SW,  i.e.,  the  number  of  grains  that  a  propagating  SW  sits  on, 


varies  slightly  and  so  does  the  SW  amplitude.  This  makes  sense  because  in  a  particulate  system  there 


would  be  time  periods  in  which  the  grains  will  get  more  compressed  and  when  they  would  be  less 


compressed  against  each  other.  Increased  compression  would  mean  high  potential  energy  whereas  smaller 
compression  would  mean  high  kinetic  energy.  This  is  the  reason  why  the  SW  height  and  spatial  extent 
varies  in  the  KE  vs.  grain  number  snapshots  in  Fig.  1.  Thus  we  focus  on  the  average  width  of  the  SW,  and 


Numerical  results  show  that  as  the  exponent  n  decreases  (but  still  greater  than  2)  the  wave  hump  tends  to 
flatten  out,  i.e.,  the  width  of  the  SW  increases  (see  Fig.  2).  In  our  numerical  studies  we  choose  to  keep  the 
total  energy  of  the  system  and  hence  of  the  SW  to  be  the  same  in  every  case  we  examine.  The  kinetic  and 
the  potential  energies  of  the  solitary  wave  then  end  up  being  related  by  the  Virial  theorem. 

The  Virial  Theorem  is  generally  valid  for  classical  energy  conserved  systems.  Virial  Theorem 
[Goldstein,  2002]  says  that 
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(4) 


<  KE  >= 


< 


If' 


r  >= 


1  y-,dV 

—  <  >  — 

2  rac.- 


X,  > 


where  F(  is  the  force  on  particle  i  located  at  rt  and  <  >  signifies  a  time  average. 

In  the  granular  alignment, 

V~S- ~(u,-uM)\ 

and  for  for  generalized  Hertz -type  potentials 

(KE)  =  (EE)  =  (5) 

77  +  2  77  +  2 

and  these  two  results  have  been  confirmed  by  all  of  our  simulations  (  see  Fig.  1*  below). 


F,„|,  KE  and  PF.  vs.  time 


n  =2.05  2.5  3.0  4.0  5.0 


1  126  251  376  501  626  751  876  1001  1126 

time  (us) 


Fig..  1  *  KE  and  PE  vs.  time  for  different  iz’s:  2.05,  2.5,  3.0,  4.0  and  5.0  from  left  to  right. 

Now  kinetic  energy  can  be  determined  as  a  function  of  n,  if  Etot  of  the  system  is  known.  Hence  the 

descriptions,  including  the  width,  of  the  SW  can  be  fully  obtained  through  the  numerical  studies  of  the 
system.  However,  as  we  shall  see,  the  relation  of  the  SW  width  to  its  kinetic  energy  is  not  trivial  to 
determine 


13 


for  different  n  values 


Fig.  2  Some  numerical  results  of  the  numerical  solutions  to  the  equations  of  motion  showing  the  velocity 
of  the  solitary  waves  for  different  magnitudes  of  n.  The  bigger  the  n  is,  the  narrower  the  solitary  wave  is. 


It  is  difficult  to  exactly  define  the  width  W  of  a  SW.  In  experiments  and  in  simulations,  the  width  is 
related  to  the  accuracy  with  which  measurement  is  possible.  In  continuum  theory  the  width  may  diverge. 
Our  calculations  show  that  the  grain-grain  distance  between  the  grains  inside  a  solitary  wave  may  differ 
by  as  many  as  six  decades.  It  is  hence  necessary  to  introduce  a  reasonable  cut-off  in  some  appropriate 
parameter  (displacement  or  velocity  say)  to  precisely  define  W. 

We  in  the  following  present  the  detailed  steps  for  calculating  the  width  W.  The  method  involves  the 
symmetry  property  of  SW.  It  says  a  SW  shape  is  line-symmetric  about  its  center,  see  Fig.  3.  The  vertical 
axis  represents  the  velocity  squared  of  each  grain  involved  in  the  SW. 
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Fig.  3.  Sketch  of  the  SW  being  replaced  by  the  isosceles  triangle.  Panel  A  illustrates  that 
Areal  +  Areal  =  Areal  where  Areal  is  very  small  compared  with  total  area  of  the  triangle.  One  can 
always  construct  such  an  isosceles  triangle  whose  area  is  equal  to  that  under  the  solitary  wave  curve. 


Virial  Theorem  (Eq.  (5))  allows  us  to  write 
Area  of  the  triangle  x  m/2  =  KE  totai  =  n  E  tota!/  (n+2),  (6) 


where  m  is  the  mass  of  each  grain. 

Expressed  now  in  2-D  coordinates  x-y,  the  area  of  the  triangle  is  (see  Fig.  3 -panel  B) 

A  =  \B(\ wvmax )  =  ^  BmvL,  >  (?) 

which  is  the  total  kinetic  energy  of  the  solitary  wave,  i.e. 


1 

4 


Bmv„ 


co  1 

=  KE=  Y,2mv‘2 


2  1 

£  mv: 
w  2 


(8) 
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where  B  =  W  +  1 .  Making  use  of  equation  (4)  and  the  definition  of  B  =  W  +  1 ,  we  have  the  width 


expression, 


W  = 


4/2 


n  + 2 


VWVmax7 


-1. 


(9) 


Eq.  (9)  allows  us  to  find  the  width  of  a  SW  for  a  given  index  n;  the  known  total  energy  of  the  system,  and 
the  center  grain’s  velocity  (which  is  the  maximum  velocity)  as  obtained  from  the  numerical  solutions. 

1  2 

Recall  that  the  total  energy  of  the  system  is  Etot  =  —  m  v0  ,  where  v0  is  the  initial  velocity  imparted  to 
the  first  grain  at  t  =  0,  therefore  one  can  also  write  the  width  equation  in  terms  of  v0  as  follows 


W  = 


2/2 


n  +  2 


(  2  A 


>  v  , 

V  max  J 


-1. 


(10) 


Eq.  (10)  provides  the  relation  between  the  SW  width  and  the  power  law  potential  exponent  n  given  the 
initial  impulse  velocity  v0  and  the  maximum  velocity  vmax  of  the  center  grain  of  the  grains  that  carry  the 
SW. 

Once  v  max  is  known  from  the  numerical  simulation  result,  the  width  of  the  SW  for  a  given  potential 
component  n  can  be  obtained.  This  process  can  be  repeated  indefinitely  until  all  studies  of  widths  for 
different  magnitudes  of  n  are  completed.  The  width  depends  upon  the  magnitude  of  exponent  n  in  the 
following  way:  n  — >  oo,  W  — »  1 ,  i.e.,  rigid  hard  potential  limit;  and/2  — »  2  ,  solitary  wave  breaks  down, 
i.e.,  harmonic  potential  limit. 

Our  exhaustive  numerical  studies  (see  Fig.  4)  of  such  physical  systems  have  shown  that  the  width 
depends  solely  upon  the  potential  exponent  n,  and  is  consistent  with  a  scaling  law, 

W  =  A0(n-2)-a  +1,  (11) 

in  which  0.3301,  and  the  pre-factor  T0=4.5262.  A0  is  a  constant  that  only  depends  on  the  cutoff 

threshold  energy.  We  studied  W  for  various  n  values.  Our  simulations  confirm  that  W  °c  (n  —  2)  0  3301  or 
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W  oc  [n  —  2)  13  where  the  index  1/3  is  approximate.  We  are  unable  to  further  clarify  whether  the  actual 
index  is  1/3  or  not. 


\n(lV-\ )  vs.  ln(«-2) 


-10  -8  -6  -4  -2  0 


ln(n-2) 

Fig.4.  We  show  a  ln-ln  plot  of  W-l  vs.  n-2.  The  width  W  depends  on  n  as:  W  =  4.5262 (n  —  2)  1/3  + 1 

Analytical  Study 

In  2001,  Sen  and  Manciu  proposed  a  series  solution  for  the  equation  of  motion.  The  solution  is  as 
follows: 

«,.(/)  =  m(z;.,/)  =  n(zi  - ct )  =  u(a)  =  A(j>n{a )  (12) 

where  a  =  z  -ct,  c  =  velocity  of  SW.  Then  we  define 

“(a)=4^»M+1l  ^(a)=_tanhf'^^]  (13) 

2  V  1  7 

and 
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(14) 


/>')=!> 

9=0 


2q+\ 


{11)2 


l2q+\ 


where  z'=  z-ct 


So  knowledge  of  c  2q+i  («)  solves  u  (z-ct)  which  reveals  the  displacement  of  every  grain  within  the  SW 
and  hence  solves  the  equation  of  motion.  u(z  —  ct )  can  hence  be  found  and  <KE>  and  <PE>  can  in 


principle  be  obtained.  The  following  table  lists  the  constants  C  1-..5  for  different  n  values  [Sen  2001].  The 


constant  CA  = 


me: 


na. 


rA^-2 


where  a  is  from  Eq.  (1). 


Table  1  Coefficients  C  values  for  different «’  s 


n=2.2 

2.35 

2.5 

3.0 

4.0 

5.0 

CO 

0.8709 

0.6908 

0.85852 

0.9445 

1.3323 

2.0517 

Cl 

1.643 

2.3171 

2.3953 

3.0168 

3.5646 

3.79001 

C3 

0.08223 

0.2364 

0.26852 

0.5971 

1.331 

2.177 

C5 

0.000326 

0.003407 

0.006134 

0.0376 

0.0676 

0.0665 

In  the  following  we  report  analytical  results  of  SW  widths  dependence  on  n  and  how  we  construct  an 
analytical  expression  of  this  SW. 

We  will  now  analytically  obtain  the  average  potential  energy  <PE>  of  the  SW.  The  displacement  of 
involved  grains  from  equilibrium  positions  has  the  form  [Sen  2001], 


ut(zi  -0 


1 


1  +  e 


fM-u)  ’ 


(15) 


where  fn  (zi  -  tf  )-Z  C,  1  (n)(zi  —ti)lq"  is  a  controlling  factor  containing  all  the  information  of  this 

?=o 

type  of  solitary  wave,  and  zt  [  =  □  □  D±D  □  DiDDD  □  □  □□  D±D  W-  1)  /  2  ]  denotes  discrete 


t 


dimensionless  coordinate  recording  the  grain  number,  /( is  a  characteristic  time,  i.e.,  t;  =  -L^~  (t0  =time 


for  solitary  wave  to  travel  one  grain  diameter)  which  is  not  necessary  to  be  discrete.  Knowledge  of  the 
coefficients  C0 ,  C1 ,  C3 ,  C5  ,•  •  •,  [6]  will  completely  solve  the  problem  of  pulse  propagation  for  any  system 
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supporting  this  type  of  SW.  One  must  resort  to  numerical  methods  for  computing  these  coefficients.  For  a 
certain  time  th  say  0  for  convenience  we  choose  ut  (z. )  to  be  centered  at  the  origin  so, 


,,(z')=TT^ 


(16) 


which  leads  to  the  overlap  of  Dy+i  of  two  adjacent  grains, 

Si,M  =Ui(Zi)-Ui+l(Zi+l)  (l7) 

Hence  the  potential  energy  is, 

w 

00  2 

PE)=  YsanSi,M  ~  Yjan$i,i+ 1"  >  (18) 


W 

2 


where  a.,  = 


f  Y  ^ 

'  R3-"  ) 

V 1  —  cr2  J 

[r-3n(n-  1)J 

,  Fx  (n  -  2,2  -  n;  n  -  l;l) .  For  the  sake  of  convenience,  we  let 


6  7  7 


Vl-cr  7 


R  '  =  1 ,  hence  a„  = 


1 


2"  3n(n- 1) 


,  Fx  ( n  -  2,2  -  n;  n  -  l;l) . 


The  kinetic  energy  of  the  SW  likewise  can  be  found  as  follows, 
u.  (t. )  = - 1  ,  (note  that f„  is  an  odd  function) 


l  +  e 


which  gives, 


(19) 


/  \  dll  I  q- 0 


Z(2^  +  l )C2?+l(«)  t 


2  q 


- f,X< ) 


(l  +  e-/„M)2 

Therefore  the  kinetic  energy  is, 


(20) 


W=Zfvl2-2  fvr 

7=— 00  ^  i=~w  n  ^ 


(21) 


where  m  is  set  to  be  unity. 
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Uj  (  relative  displacement ) 


<5u+i  (overlap) 


Fig.5  The  figure  is  showing  the  relative  displacement  of  each  bead  in  a  solitary  wave,  and  the 
compressions  (overlaps)  of  pairs  of  neighboring  beads  which  constitute  a  compression  solitary  wave. 

With  the  knowledge  of  coefficients  Cl,C3,  C5 ■  •,  to  desired  accuracy,  the  SWs  potential  and  kinetic 
energies  can  thus  be  calculated  for  all  various  n  values.  Hence  the  widths  have  been  determined.  The 
accuracy  in  our  studies  is  to  one  part  in  1020.  We  calculated  8  data  sets  of  8 i  i+1  at  different  times  to  get  < 

8 ]  (+1  >  shown  in  Table  2  below. 
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Table  2.  Relative  displacements  for  different  n  values 


(au) 

n  =2.2 

n  =2.35 

n  =2.5 

n  =3.0 

n  =4.0 

n  =5.0 

z 

D □□□□□ 

0 

3.2449E-23 

7.88E-29 

1.01  IE-90 

1.8E-172 

2.167E-217 

-6 

n □□□□□ 

5.18961E-06 

7.7427E-13 

4.44E-15 

2.764E-39 

5.62E-74 

2.1527E-97 

-5 

D □□□□□ 

0.000719994 

7.0716E-07 

1.21E-07 

1.26E-15 

4.12E-28 

3.2861E-38 

-4 

n □□□□□ 

0.018086981 

0.0013119 

0.000796 

6.061E-06 

2.19E-09 

1.6602E-12 

-3 

Q □□□□□ 

0.132344218 

0.07065122 

0.063974 

0.0252896 

0.006942 

0.00239133 

-2 

n  □□□□ 

0.348843614 

0.42803616 

0.43523 

0.4747043 

0.493058 

0.49760867 

-1 

Q  □□□ 

0.348843614 

0.42803616 

0.43523 

0.4747043 

0.493058 

0.49760867 

0 

n  □□□ 

0.132344262 

0.07065122 

0.063974 

0.0252896 

0.006942 

0.00239133 

1 

Q  □□□ 

0.018086937 

0.0013119 

0.000796 

6.061E-06 

2.19E-09 

1.6603E-12 

2 

n  □□□ 

0.000719995 

7.0716E-07 

1.21E-07 

1.258E-15 

4.12E-28 

3.2861E-38 

3 

Q  □□□ 

5.18942E-06 

7.7434E-13 

4.44E-15 

2.764E-39 

5.62E-74 

2.1527E-97 

4 

n  □□□ 

0 

3.2449E-23 

7.88E-29 

1.01  IE-90 

1.8E-172 

2.167E-217 

5 

width 

9 

7.8 

7 

5.8 

4.6 

4.4 

To  ensure  the  width  found  by  the  guidance  of  cutoff  ratio  of  8tail  /  8 max  <10  8  is  correct,  W  must  also 
satisfy  the  following  relations 


<KE>  = 

w 


f  n  ^ 

n  +  2 


E,ot,  and  <  PE  >  w  = 


r  2  " 

n  +  2 


E,o,  or 


<  KE  >  n 


<  PE  > 


(22) 


The  following  table  lists  the  KE,  PE  and  Etot  for  different  n  values. 


Table  3.  KE  and  PE  for  different  n  values 


n 

2.2 

2.35 

2.5 

3 

4 

5 

<KE>  (arb.u) 

0.29131754 

3 

0.52372503 

0.96139851 

8 

1.34954430 

3 

2.09294106 

5 

2.42609335 

9 

expected 

<KE>/<PE> 

1.1 

1.175 

1.25 

1.5 

2 

2.5 

<PE>  (arb.u ) 

0.26546468 

7 

0.44678467 

0.77095005 

0.90083918 

4 

1.04896212 

9 

0.97274790 

9 

calculated 

<KE>/<PE> 

1.09738717 

3 

1.17220903 

1.24703087 

9 

1.49809680 

4 

1.99524940 

6 

2.49406175 

8 

Etot  (arb.u) 

0.55678223 

0.9705097 

1.73234856 

8 

2.25038348 

7 

3.14190319 

5 

3.39884126 

7 

Prefactor(  arb.u 
) 

0.20872971 

4 

0.22953622 

0.528041 

4.20465836 

2 

0.47439964 

9 

0.51389658 

6 
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It  is  clearly  seen  that  the  Eqs.  (22)  are  satisfied  closely.  The  plot  of  expected  <KE>/<PE>  vs.  calculated 


<KE>/<PE>  shows  the  same  agreement  (see  Fig.  6) 


<KE>/<PE>  vs.  n 


n 


Fig.  6.  Both  calculated  and  expected  <KE>/<PE>  match  closely  vs.  different  n’s. 

It  turns  out  that  (see  the  plot  of  W vs.  n  shown  in  Fig.. 7)  dependence  of  W on  n  is 

W  =  4.745(n  -  2r°'3307  +  1  ,  (22) 

which  agrees  with  the  results  obtained  by  the  numerical  approach. 

The  solution  is  approximate  because  Cl,C3,C5,---  used  is  an  infinite  series. 
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ln(W-l)  vs.  ln(/i-2)  from  analytical  results 


Fig.7.  The  widths  of  PE  SW  being  calculated  and  plotted  against  n.  The  ln-ln  plot  of  W- 1  vs.  n-2  shows 
the  width:  W  =  4.745(«  -  2)~0J307  +  1  ,  which  provides  a  good  agreement  with  the  scaling  law. 
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Fig.8.  Analytical  results  of  SWs  (potential  type)  for  set  of  different  potential  exponents. 
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grain  number 


Fig.9.  Log  scale  potential  type  SWs  for  set  of  different  potential  exponents. 

Summary 

To  summarize  this  work,  we  conclude  that  the  width  of  the  SW  in  granular  media  exhibits  an  intrinsic 
power  law  dependence  on  the  potential  exponent  n,  i.e.,  W  qc  (n  -  2)  13 +1.  We  have  obtained 
agreement  between  the  numerically  and  theoretically  obtained  values  of  this  exponent.  In  monodispersed 
chains,  n  is  the  only  parameter  that  controls  the  width  of  the  SW.  The  physical  meaning  of  the  scaling 
law  indicates  that  the  width  of  the  SW  decreases  much  faster  than  the  increase  of  the  exponent  n  which  is 
solely  due  to  the  nature  of  power  law  form  of  the  potential.  Therefore,  when  n  increases,  the  energy 
transfer  from  one  grain  to  the  next  becomes  a  progressively  faster  process.  Such  energy  transfer  tends  to 
happen  as  a  “one-shot  deal”.  The  exponent  in  scaling  law  (-1/3)  must  hence  be  greater  then  -  1  and  less 
than  zero. 
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Task  (iii):  In  recent  work  using  compressed  granular  alignments  using  reflective  walls  we  have  shown 
that  there  exists  a  critical  loading  where  solitary  waves  and  acoustic  waves  can  coexist  with  the  kinetic 
energy  of  the  solitary  wave  being  slightly  larger  than  the  average  kinetic  energy  of  the  acoustic  waves. 
Here  we  will  attempt  to  develop  a  system  where  the  solitary  wave  can  be  submerged  under  the  noise  by 
suitably  modifying  the  system  Hamiltonian  in  such  a  way  that  the  added  noise  will  not  modify  the 
coexisting  phase  ( Submerged  Solitary  Wave  Problem). 


We  came  to  understand  that  the  real  challenge  was  to  make  a  stable  solitary  wave  within  the 
quasi-equilibrium  phase  and  this  is  what  is  accomplished  below.  This  work  extends  earlier  work  which 
hinted  at  a  stable  solitary  wave  in  the  quasi-equilibrium  phase  at  critical  loading  where  the  harmonic  and 
nonlinear  forces  are  equally  strong  competitors.  The  work  suggests  that  at  such  critical  strengths, 
combination  of  nonlinear  and  linear  forces  can  lead  to  unexpected  nonlinear  structures.  The  solitary  wave 
can  be  submerged  within  noise  without  decay  but  this  is  not  possible  to  do  by  using  a  single  Hamiltonian. 
We  present  below  in  details  our  work  on  these  stable  solitary  waves  in  the  highly  noisy  quasi-equilibrium 
phase. 

Long  lived  solitary  wave  in  a  ID  granular  chain  (Yoichi  Takato,  Edgar  Avalos,  Surajit  Sen,  SUNY 
Buffalo) 

Abstract  -  Kinetic  energy  uctuations  of  a  non-dissipative  ID  granular  chain  held  between 
perfectly  reecting  walls  with  various  precompressions  are  numerically  investigated.  The  dynamics  of  the 
system  is  explored  for  three  cases,  (i)  for  a  weakly  precompressed  chain  which  admits  only  solitary  waves 
that  break  down  into  secondary  waves  under  wall  collisions,  which  is  consistent  with  other  ndings  in  the 
literature,  (ii)  for  a  strongly  precompressed  chain  which  exhibits  mostly  acoustic  waves,  and  (iii)  for  a 
chain  with  intermediate  precompression.  Case  (iii)  accommodates  a  nearly  stable  solitary  wave  that 
travels  unaffected  through  the  acoustic  waves  and  bounces  back  and  forth  between  the  walls  for 
extremely  long  times. 

Nesterenko  was  the  first  to  report  the  existence  of  a  traveling  solitary  wave  (SW),  i.e.,  a  non-dispersive 
bundle  of  energy,  in  a  ID  granular  chain  [1].  Since  then,  many  theoretical  and  experimental  studies  on 
SWs  have  been  carried  out  [2-4].  These  studies  have  shed  light  on  the  intriguing  characteristics  of  SWs  in 
granular  chains  where  grains  repel  via  the  strongly  nonlinear  Hertz  potential  [5].  One  of  the  properties  of 
SWs  that  was  shown  by  computational  simulations  is  as  follows:  a  SW  created  out  of  _  function 
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perturbation  to  an  edge  grain  of  a  chain  breaks  down  into  small  secondary  SWs  due  to  interactions  with 
walls  or  other  SWs  [6],  Although  the  secondary  SWs  carry  much  lower  energy  compared  to  the  initial 
SW,  they  were  experimentally  found  by  Job  et  al.  [4],  This  breakdown  of  SWs  eventually  leads  the  chain 
to  an  interesting  state  as  we  shall  see  below. 

Suppose  that  only  one  SW  is  created  by  a  single  initial  impulse  in  the  system.  Now  imagine  that 
through  wall  collisions  and  collisions  with  one  another  many  secondary  SWs  are  born  out  of  the  single 
original  SW.  One  would  expect  that  eventually  the  chain  will  reach  a  state  that  will  be  full  of  secondary  or 
low  energy  SWs.  In  this  state,  the  secondary  SWs  will  continue  to  exchange  energy  with  one  another  via 
collisions.  Energy  conservation  demands  that  eventually  the  secondary  SWs  will  break  down  and  reform 
at  the  same  rate  in  the  system.  One  may  think  that  this  is  an  equilibrium  state  by  analogy  with  the 
equilibrium  state  of  a  spring-mass  system  whose  potential  is  described  by  a  harmonic  potential.  However, 
it  turns  out  that  the  actual  situation  is  quite  different.  For  instance,  the  velocity  distribution  in  this  state  is 
a  Gaussian,  whereas  the  equipartition  theorem  does  not  hold  [7],  Further,  in  some  cases  the  system's 
evolution  appears  to  have  no  memory  of  the  initial  conditions  whereas  in  others  they  do.  In  addition, 
when  one  investigates  the  fluctuations  in  kinetic  energy  for  the  Hertzian  system  in  the  equilibrium-like 
state  one  finds  that  the  fluctuations  are  larger  than  those  for  a  harmonic  system  [7].  The  conclusion  -  this 
state  is  not  the  same  as  the  equilibrium  state  in  a  harmonic  chain.  Therefore,  this  state  was  named  as  the 
quasi-equilibrium  state  to  distinguish  it  from  the  ordinary  equilibrium  state  [7,  8]. 

The  dynamics  of  anharmonic  systems  has  fundamental  differences  with  respect  to  those  of 
harmonic  systems  [9].  In  particular,  the  Hertzian  system  has  an  important  property  which  comes  from  the 
discreteness  of  the  system;  the  weakly  compressed  system  does  not  allow  acoustic  waves  to  propagate, 
these  systems  have  been  referred  to  as  being  in  sonic  vacuum  [1],  On  the  contrary,  the  precompressed 
chain  can  admit  acoustic  waves  [  1 0]  and  strong  precompression  may  suppress  the  nonlinearity  and  lead  to 
a  state  with  equipartitioned  energy  in  the  system  like  in  a  typical  harmonic  system.  The  role  of 
intermediate  precompression  in  the  Hertzian  system  is  not  well  understood. 

In  this  paper  we  present  our  numerical  study  of  kinetic  energy  fluctuations  of  a  finite 
nondissipative  linear  chain  made  of  Hertzian  grains  with  precompression.  The  precompression  ranges 
from  weak  enough  so  as  to  achieve  sonic  vacuum  to  strong  enough  so  as  to  have  acoustic  waves  in  the 
system.  We  present  evidence  to  show  the  existence  of  a  striking  long-lived  SW  at  intermediate 
precompression. 

Model  and  Method:  We  consider  a  ID  chain  composed  of  N monodispersed  spherical  grains  of 
radius  R  with  rigid  walls  at  both  ends  of  the  chain.  These  grains  are  barely  in  contact  with  one  another  at 
zero  precompression.  When  two  adjacent  spherical  grains  are  compressed,  the  repulsive  potential  is 
described  by  the  Hertz  potential  [5], 
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(1) 


K<w)  =  { 


a^i,i+ 1  0^i,i+ 1  ^  ® 
0  ( otherwise ) 


where  n  =  5/2  and  <5(  i+1  =  2R  —  (iq+1  —  ut)  is  the  overlap  between  grain  i  and  i  +  1  and  iq,  ui+1 
describe  the  positions  of  the  centers  of  grain  i  and  i  +  1,  respectively.  The  prefactor  a  is  given  by 
a  =  (2/SD)^J R /2  ,  where  D  is  determined  by  the  grain’s  material  properties;  Poisson  ratio  a  and 
Young’s  modulus  Y  give  D  =  3(1  —  o2)/2Y.  If  the  separation  of  two  neighboring  grains  is  greater  than 
2 R,  they  do  not  have  physical  contact  and  hence  no  force  acts  on  the  grains.  For  this  study  we  set 
a  =  4.136  X  107N/m3  2  and  R  =  0.5mm.  This  value  of  a  corresponds  to  materials  such  as  glasses. 

The  Hamiltonian  is  described  by 

=  Z?=i||  +  Zf=iK<W)>  (2) 


where  pi  is  the  momentum  of  grain  i,  m  is  the  mass  of  a  grain  and  V (Si  i+1)  is  the  potential  energy  fue  to 
the  overlap.  The  equation  of  motion  of  grain  i  which  is  not  next  to  a  wall  is  given  by 

m^T  =  nai\-A  +  ui- 1  -  “;]n-1  -  [A  +  ut  -  U(+1]n_1  (3), 
where  the  symbols  are  as  defined  above  and  A  is  the  precompression. 

A  velocity  perturbation  is  given  to  one  end  of  the  chain  at  t  =  0  and  it  eventually  develops  into  a 
traveling  SW.  The  impact  velocity  is  set  to  be  vimp  =  9.899  X  10-2  m/s  for  all  cases  studied  here.  We 
solve  the  equation  of  motion  using  the  velocity  Verlet  algorithm.  The  integration  time  step  is  At  =  10~9s. 
Now,  we  define  a  time  averaged  fluctuation  of  the  kinetic  energy  T  of  the  chain  as  [7] 


<F(t)>=Sa""|T,(t/r>'  U) 

where  Tt{t)  is  the  instantaneous  kinetic  energy  of  the  chain,  <  T  >  is  the  average  kinetic  energy  at  t  -*  oo, 
,Nm  is  the  total  number  of  time  steps,  and  Ft(t)  is  the  instantaneous  fluctuation  of  the  kinetic  energy. 

Fluctuations  of  kinetic  energy  of  a  granular  chain:  We  investigate  how  kinetic  energy  fluctuations 
behave  as  the  precompression  applied  to  the  chain  is  varied  for  three  different  system  sizes  N-  50;  100, 
and  500.Fig.  1  shows  the  kinetic  energy  fluctuations  obtained  by  our  simulations.  In  each  system  size,  the 
fluctuation  is  normalized  by  the  value  of  the  fluctuation  at  no  precompression.  In  addition,  the 
precompression  A  is  nondimensionalized  by  grain  radius  R  =  0:5  mm. 
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There  are  two  distinct  states  at  high  and  low  precompressions  in  _g.  1.  At  the  highest  and  lowest 
precompressions  the  fluctuations  are  almost  at,  however,  the  levels  of  the  fluctuations  are  different  in  the 
two  cases.  The  fluctuation  for  low  precompression  is  20%  larger  than  that  for  high  precompression  and 
hence  shows  a  clear  difference  between  the  two  states.  This  is  consistent  with  the  results  reported  by  Sen 
[7]  that  quasi-equilibrium  state  achieved  in  a  weakly  precompressed  chain  in  small  systems  has  more 
fluctuations  by  12  {27%  compared  to  the  equilibrium  state.  The  equilibrium  state  corresponds  to  the  state 
established  in  the  strongly  precompressed  chain  which  allows  mostly  acoustic  wave  propagation. 
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Fill.  2.  Tin-  exponent  a  for  tin-  A  MB  syitim 

Another  common  feature  to  all  the  panels  in  fig.  1  is  that  the  uctuations  possess  a  broad  peak  at  an 
intermediate  precompression  regime  between  1 0~4  and  1 0‘2.  Our  simulation  suggests  that  this  broad  peak 
in  the  fluctuations  has  to  do  with  the  fact  that  the  fluctuations  are  yet  to  reach  a  steady  state  during  the 
calculation  of  the  fluctuations.  Further,  the  inability  to  reach  a  steady  state  is  because  the  kinetic  energy  at 
the  precompression  of  interest  changes  extremely  slowly  in  time.  This  is  confirmed  in  Fig.  2,  which 
shows  an  exponent  a  such  that  max  {7)(t)}  oc  t~a .  The  maximum  kinetic  energy  of  the  weakly 
compressed  chain,  which  should  be  attributed  to  the  initial  solitary  wave,  decays  very  quickly,  and  its 
decay  rate  does  not  depend  on  the  precompression,  as  shown  in  Fig.  2  as  a  large  value  that  remains  fixed 
for  A<  2  X  10-4.  By  contrast,  the  maximum  kinetic  energy  of  the  chain  with  intermediate 
precompression  decays  slowly,  and  the  rate  becomes  smaller  with  the  increasing  precompression 
around2  x  10~4  <  A<  1.5  x  10-2.  (see  the  negative  slope  in  fig.  2).  It  indicates  that  uctuations  of  the 
kinetic  energy  at  intermediate  precompression  could  possibly  be  altered  if  the  runs  were  even  longer. 
Otherwise  the  fluctuation  for  each  system  size  may  converge  to  the  smooth  solid  line  shown  in  fig.  1 , 
which  is  obtained  using  points  not  belonging  to  the  peak.  To  run  until  convergence  of  the  fluctuations, 
however,  costs  too  much  time  as  the  exponent  values  show.  For  A=  9  X  10-4,  we  ran  the  simulations 
about  thirty  times  longer  (a  month  in  real  time)  than  those  for  weakly  compressed  cases.  Finally,  for  the 
strongly  compressed  chain,  the  exponent  -»  oo  because  an  initial  solitary  wave  which  causes  large 
uctuations  may  not  exist.  The  system  reaches  a  steady  state  at  the  very  beginning  of  the  runs,  which  will 
be  presented  below  in  fig.  3(e),  and  the  average  maximum  kinetic  energy  remains  unchanged  in  time. 

So  far,  we  have  looked  at  the  kinetic  energy  fluctuations  of  the  system  in  the  presence  of 
precompression.  The  results,  large  fluctuations  are  seen  for  low  precompression  and  small  fluctuations  are 
found  for  high  precompression.  These  findings  are  not  surprising  because  of  the  properties  of  the  solitary 
waves  in  question.  These  waves  break  down  into  smaller  waves  and  eventually  form  a  steady  state,  which 
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leads  to  the  quasi-equilibrium  state  for  low  precompression  and  at  large  precompressions  the  acoustic 
waves  dominate,  and  these  lead  to  the  system  being  in  a  state  with  equipartitioned  energy.  What  is 
unexpected  is  the  odd  slowing  down  of  the  system  towards  a  steady  state.  This  form  of  behavior  is 
reminiscent  of  critical  slowing  down  in  systems  undergoing  phase  transition  [11]  and  we  focus  on  this 
intriguing  behavior  below. 
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A  space-time  plot  of  kinetic  energy  shown  in  gray  scale  (a  “zigzag"  plot)  helps  us  easily  visualize 
traveling  SWs  or  acoustic  waves  in  the  chain  in  order  to  explore  where  the  decay  rate  difference  comes 
from.  The  darker  spots  in  the  plot  represent  higher  kinetic  energy.  Such  a  spot  may  display  a  bundle  of 
kinetic  energy  moving  through  the  chain  or  how  kinetic  energy  spreads  in  the  chain.  A  SW,  a  bundle  of 
kinetic  energy,  is  shown  as  a  sharp  line  in  the  zigzag  plot.  If  a  SW  propagates  from  one  end  of  the  chain 
to  the  other  end,  the  plot  shows  a  straight  line  with  a  certain  slope,  which  is  the  velocity  of  the  SW.  If 
kinetic  energy  is  equally  distributed,  i.e.,  an  equilibrium  state  is  established,  the  plot  shows  a  uniform 
gray  pattern.  The  time  shown  in  fig.  3  is  normalized  by  the  period  of  a  harmonic  oscillator  tharm  = 

2n^m/k,  where  k  is  the  spring  constant  obtained  by  letting  k  =  a  with  Si  i+1  =  R. 


32 


Upon  investigating  the  zigzag  plots  for  the  entire  precompression  window  for  an  N=  500  system,  we  find 
that  the  patterns  in  the  plots  can  be  classified  into  three  distinct  types  as  shown  in  fig.  3.  Fig.  3(a)  and  (b) 
represent  precompression  ofA<  2  X  10-5  (corresponding  to  weak  precompression),  fig.  3(c)  and  (d) 
represent  2  X  10-4  <  A<  10-2  (corresponding  to  intermediate  precompression),  and  fig.  3(e)  and  (f) 
represent  A>  10-2  (corresponding  to  strong  precompression).  A  zigzag  plot  for  the  weakly  compressed 
chain  in  fig.  3(a)  shows  that  a  clear  initial  SW  depicted  by  a  thin  solid  line  at  the  beginning  (t  <  300) 
travels  in  the  chain;  it  starts  to  decay  shortly  after  t  =  300;  it  decays  quickly,  breaking  down  to  secondary 
SWs  around  t  =  500.  In  the  case  of  the  strongly  compressed  chain,  the  zigzag  plot  is  drawn  in  fig.  3(e).  It 
exhibits  the  behavior  that  is  characteristic  of  that  of  sound  waves.  The  initial  perturbation  spreads  into  the 
chain,  and  more  than  half  of  the  grains  in  the  chain  have  some  kinetic  energy  by  the  time  the  wave  front 
reaches  the  opposite  wall.  Therefore,  it  cannot  be  a  SW.  Moreover,  at  a  later  time  the  plot  pattern  turns 
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uniformly  gray,  meaning  that  the  system  is  now  in  a  state  with  equipartitioned  energy  with  seemingly  no 
dependence  on  initial  conditions  and  shows  a  Gaussian  distribution  of  velocities  centered  around  zero. 
Lastly,  for  the  chain  with  intermediate  precompression,  a  single  SW  and  harmonic  waves  coexist  as  can 
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be  seen  in  _g.  3(c).  A  SW  represented  by  a  sharp  line  in  the  plot  is  observed  in  fig.  3(a),  and  the  gray 
pattern  similar  to  the  one  shown  in  fig.  3(e)  is  also  observed.  Surprisingly,  the  SW  lives  for  an  extremely 
long  time,  until  at  least  t  =  6.9  X  104,  which  is  the  end  of  the  run.  The  SW  does  not  attenuate  much  even 
though  it  repeatedly  collides  with  the  walls.  The  precompression  exerted  on  the  system  seems  to  play  an 
important  role  in  suppressing  the  mechanisms  that  are  associated  with  the  decay  of  the  SW  via 
secondarySW  formation  and  this  seems  to  lead  to  the  longevity  of  the  SW. 

Static  and  dynamic  potential  energy,  and  crossover  precompression  Ac:  The  precompression  of 
the  system  is  the  only  parameter  varied  for  each  system  size  in  our  study  of  fluctuations  of  kinetic  energy. 
Impact  velocity,  grain  size,  grain  mass,  and  all  the  material  properties  associated  with  this  investigation 
are  kept  the  same.  The  precompression  is  imparted  to  the  system  as  static  potential  energy.  Hence,  the 
effect  of  precompression  is  similar  to  that  of  an  external  field  that  affects  the  system's  dynamics.  As  we 
have  seen  above,  the  SW-acoustic  wave  coexistence  state  emerges  at  a  certain  precompression  value  and 
this  may  be  a  special  state  of  the  system.  As  argued  extensively  elsewhere  [8],  at  zero  precompression 
approximately  5/9  of  the  system’s  energy  turns  out  to  be  kinetic  and  the  remaining  4/9  to  be  potential. 
Such  a  system  when  subjected  to  a  single  velocity  perturbation  forms  a  solitary  wave  pulse  that  carries  all 
of  the  system's  energy.  Precompression  changes  this  behavior.  Upon  precompression,  the  system  can 
possess  more  potential  energy  and  when  the  amount  of  potential  energy  due  to  loading  equals  that  due  to 
the  Hertzian  interaction,  a  new  state  forms.  In  this  state,  the  SW  is  present  but  the  extra  potential  energy  is 
used  to  generate  acoustic  waves  which  coexist  alongside  the  SW.  In  _g.  4(a)  the  potential  energy  and  the 
kinetic  energy  of  the  chain  are  shown  as  a  function  of  the  precompression.  The  potential  energy  at  high 
precompression,  A>  10~3,  A>  10-3,  is  dominated  by  the  precompression  and  that  at  low 
precompression,  A<  10-4,  is  dominated  by  the  Hertzian  interacton.  The  coexistence  regime  begins 
where  the  the  potential  energy  due  to  precompression  becomes  comparable  with  that  due  to  Hertzian 
interaction.  We  assume  that  a  way  to  estimate  the  strength  of  loading  where  the  SW  and  the  acoustic 
waves  coexist  is  by  defining  a  Acin  such  a  way  that  potential  energy  from  precompression  equals  the 
potential  energy  from  Hertzian  interactions.  The  details  of  this  derivation  are  shown  in  Appendix  A.  The 
result  of  the  crossover  point  is  given  by 
/  2  \l/n 

*c={^<T>‘)  •  <5> 

where  <  T  >t  is  the  time  averaged  kinetic  energy  of  the  chain.  We  perform  a  numerical  calculation  of  the 
crossover  point  using  the  equation  and  the  values  obtained  in  our  simulation  for  N  =  500,  <  T  >t= 

5.2  X  10_1  from  our  simulations,  then  we  get  Ac=  1.8  X  10-4.  It  agrees  with  our  simulation  result  of 
1.9  X  10-4.  Furthermore,  the  peak  shifts  shown  in  fig.  1  can  be  explained  by  eq.  (5)  because  the  system 
size  N  is  increased,  then  the  crossover  point  becomes  small,  which  also  explains  our  simulation  results. 
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The  acoustic  wave  velocity  goes  up  with  increasing  precompression  in  the  coexistence  regime  as 
shown  in  fig.  4(b)  and  it  approaches  the  solitary  wave  velocity  almost  independent  of  precompression. 
The  solitary  lives  longer  and  longer  as  the  two  velocities  come  close  together.  Our  investigations  suggest 
that  the  wavelength  of  the  acoustic  waves  is  slightly  larger  than  or  of  the  same  length  as  the  typical  width 
of  the  SW.  We  conjecture  that  the  SWs  are  hence  “invisible"  to  the  acoustic  waves  and  this  may  be  why 
the  acoustic  waves  cannot  break  down  the  SWs.  However,  when  the  precompression  is  further  increased, 
the  situation  changes.  The  wavelength  of  the  acoustic  waves  become  smaller  and  these  waves  can  now 
\pry  open"  the  SWs  which  are  quite  a  bit  larger  than  them,  thereby  making  them  unstable.  This  is  when 
the  system  shows  acoustic  like  behavior. 

Discussion 

In  this  study,  we  have  explored  the  problem  of  impulse  propagation  in  precompressed  grains 
that  are  held  between  rigid  walls  in  a  granular  chain.  At  vanishingly  small  and  zero  loading  we  find  that 
the  SW  breaks  down  into  the  quasi-equilibrium  state  that  is  characteristic  of  strongly  nonlinear  systems  as 
discussed  extensively  in  the  literature  now.  It  may  be  noted  that  in  quasi-equilibrium,  the  energy  is  not 
equipartitioned  in  the  system.  The  velocity  distribution  is  Gaussian  and  in  the  case  of  the  single  impulse 
problem  there  appears  to  be  no  dependence  on  initial  conditions  [7].  When  the  precompression  is 
sufficiently  large  such  that  the  grains  can  only  oscillate  harmonically  about  their  squeezed  positions,  the 
system  behaves  like  a  harmonic  chain  with  N  modes  where  N  is  the  number  of  particles.  Such  a  system  of 
course  satisfies  the  equipartition  theorem,  shows  a  Gaussian  distribution  of  velocities  and  also  has  no 
dependence  on  initial  conditions  and  can  be  regarded  as  in  an  equilibrium  state.  In  the  N  -*  oo  limit, 
though  the  distinction  between  the  quasi-equilibrium  and  equilibrium  phases  becomes  vanishingly  small. 

The  striking  result  of  this  study  is  how  the  SW  and  acoustic  waves  coexist  in  this  system  for 
values  of  precompression  between  the  two  extremes.  To  our  knowledge  such  a  coexistence  phase  has  not 
been  seen  in  granular  chains  in  earlier  work.  Here  we  find  that  there  is  a  certain  value  of  precompression, 
that  depends  on  the  Hertz  potential  and  the  system  size,  around  which  both  the  SW  and  acoustic  waves 
coexist.  This  crossover  point  appears  to  be  when  the  average  grain-grain  separation  distance  and  the 
precompression  become  equal.  Naturally  in  such  a  limit,  the  equation  of  motion  becomes  identical  to  that 
of  the  equation  of  motion  at  zero  precompression  upto  a  coupling  constant.  The  wavelength  of  the 
acoustic  waves  seems  comparable  to  the  width  of  the  SW.  However,  it  remains  to  be  understood  why 
these  SWs  do  not  break  down  significantly  via  secondary  SW  formation. 

The  research  has  been  supported  by  the  Army  Research  Office. 
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Task  (iv):  Here  we  will  focus  on  studies  where  we  shall  explore  heterogeneous  granular  alignments  with 
the  capability  of  creating  significant  hot  spots  or  energy  collection  regions  (Granular  Hot  Spot  Creation 
Problem). 
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Figure  1:  Here  we  show  the  time  evolution  of  a  119  grain  chain  for  n=2.5  in  (a)  and  for  n=2.1  in  (b). 
Multiple  perturbations  are  initiated  at  t=0.  The  dark  patterns  that  emerge  temporarily  are  the  hot  spots 
we  were  seeking.  Further  work  to  classify  the  intensity  of  the  hotspots  is  now  under  way  for  driven, 
dissipative  systems. 


The  granular  hot  spot  creation  problem  is  very  much  in  progress  and  is  likely  to  be  for  quite  some  years. 
This  is  a  challenging  problem  that  is  similar  to  that  of  the  emergence  of  rogue  waves  in  oceans  which 
form  unexpectedly  and  are  known  to  be  a  threat  to  vessels. 

The  PI  and  his  team  acknowledges  STIR  support  for  the  work  done  so  far  which  includes 
extensive  simulations  on  the  evolution  of  multiple  perturbations  initiated  at  the  initial  instant  in  a  granular 
chain  held  between  rigid  walls.  Studies  have  been  done  for  systems  with  no  dissipation.  We  found  that 
when  the  initial  perturbations  are  initiated  in  multiple  grains,  there  is  a  strong  chance  that  regions  with 
very  high  kinetic  energy  will  emerge  unpredictably  at  late  times.  These  high  energy  regions  in  space  are 
temporary  (see  Fig.  1  above).  Thus  far  we  have  found  no  way  to  predict  the  emergence  of  such  hotspots 
though  we  find  that  many  perturbations  and  specific  boundaries  are  needed  for  these  to  likely  form. 
Substantial  progress  on  this  problem  is  described  in  the  manuscript  included  with  this  report  which 
appeared  in  Physical  Review  E  late  in  2011  (E  Avalos  et  al,  Phys.  Rev.  84,  046610  (2011)).  The  next 
stage  work  will  likely  involve  an  experiment-theory  collaboration  which  is  needed  for  a  breakthrough  in 
our  understanding  -  this  phase  is  now  under  development. 
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Task  (V):  We  will  build  on  previous  studies  on  mechanical  energy  propagation  in  confined  3D  granular 
beds  of  mono-sized  Hertz  spheres  to  extend  our  code  to  incorporate  the  presence  of  air  in  interstitial  space 

{3D  Granular  Bed  with  Air  Problem). 


In  the  STIR  proposal,  the  PI  had  claimed  that  work  will  be  started  in  coding  the  case  of  3D  granular 
systems  with  interstitial  fluids  in  them.  This  work  breaks  into  new  territory  in  modeling  granular  beds 
with  a  molecular  fluid  in  the  interstitial  spaces  where  the  calculations  are  done  at  a  length  scale  of  10s  of 
microns  to  a  few  millimeters.  Matt  Westley  (PhD  student)  and  the  PI  initially  started  this  work  in  close 
dialog  with  the  Livermore  geophysics  group’s  Dr  Otis  Walton  and  Dr  Eric  Herbold.  It  turned  out, 
however,  that  essentially  no  progress  had  been  made  on  this  problem  that  is  very  useful  beyond  the  work 
of  Professor  A.J.C.  Ladd  of  the  University  of  Florida,  which  used  hydrodynamic  assumptions  that  mimic 
a  dilute,  incompressible  gas. 

Our  work  has  hence  focused  on  the  dynamics  of  a  gas  in  the  interstitial  space  between  grains.  The  focus 
has  been  on  developing  a  Molecular  Dynamics  based  calculation  scheme  that  describes  a  Lennard-Jones 
fluid.  The  molecules  suffer  elastic  collisions  with  the  walls,  which  are  made  of  surfaces  of  the  grains  that 
define  the  interstitial  space  as  detailed  below. 

1.  Introduction 

The  properties  of  energy  transmission  through  granular  media  are  generally  not  well  understood.  We 
intend  to  build  a  computational  model  for  the  propagation  of  mechanical  energy  through  such  materials  as 
dirt,  mud,  and  sand.  This  is  to  be  accomplished  by  first  exploring  the  simplest  case  of  two  grains 
impacting  while  cushioned  by  an  interstitial  fluid  (whether  air,  water,  or  something  else)  in  full  detail.  By 
constructing  a  molecular  dynamics  simulation  from  the  ground  up,  we  should  be  able  to  extract  enough 
relevant  information  to  make  useful  approximations  for  many  grains.  These  approximations,  should  they 
be  computationally  inexpensive  enough,  will  allow  us  to  create  a  full  model  for  energy  propagation 
through  large  beds  of  grains  with  interstitial  fluids  such  as  air,  water  etc. 

2  Current  progress 

Thus  far  we  have  considered  the  first  part  of  this  study.  The  molecular  dynamics  program  we  use  is  the 
Large-scale  Atomic/Molecular  Massively  Parallel  Simulator  (LAMMPS),  developed  by 
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Sandia  National  Laboratories  [1],  With  LAMMPS,  we  have  constructed  the  grain  collision  as  follows: 
Take  a  large  collection  of  Lennard-Jones  (LJ)  particles  in  equilibrium  at  a  known  temperature,  and  cut 
two  spheres  from  this  material.  These  are  to  be  our  grains.  Then,  place  these  grains  in  a  tube  or  cylinder 
which  is  a  tight  fit  around  them.  (We  do  this  to  follow  more  closely  the  experimental  study  by  Job  et.  al. 
[2].)  Then,  we  place  an  interstitial  fluid  consisting  of  single  high-temperature  light  LJ  particles  inside  of 
this  cylinder.  The  grains  are  then  each  given  an  initial  velocity  towards  the  center  of  the  cylinder  where 
they  collide  in  a  head-on  collision. 

We  record  all  relevant  data  from  this  collision,  in  particular  the  force  between  the  two  grains  (as  measured 
by  summing  the  forces  between  each  particle  within  the  clusters)  and  the  overlap  between  the  grains  for 
many  different  runs.  Currently  a  single  run  involves  the  collision  of  two  1500-particle  grains  at  a  given 
initial  velocity,  and  takes  approximately  two  minutes  of  computer  time  to  run  with  acceptable  numerical 
error  bounds.  The  next  step  is  to  take  all  of  the  data  that  has  been  generated  in  this  way  and  analyze  it  to 
extract  the  force  law  between  the  two  spheres.  Ordinarily,  we  would  expect  a  Hertz  law  for  the  sphere 
contact,  but  the  interstitial  fluid  has  proven  to  have  a  large  effect  on  the  force  law  [2].  We  hope  to  recover 
the  linear  approximation  made  by  Job  et.  al.  [2]  for  this  collision. 
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